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Post-translational modifications of histone proteins are an important factor in epigenetic control 
that serve to regulate transcription, depending on the particular modification states of the histone 
proteins. We study the stochastic dynamics of histone protein states, taking into account a feedback 
mechanism where modified nucleosomes recruit enzymes that diffuse to adjacent nucleosomes. We 
map the system onto a quantum spin system whose dynamics is generated by a non-Hermitian 
Hamiltonian. Making an ansatz for the solution as a tensor product state leads to nonlinear partial 
differential equations that describe the dynamics of the system. Multiple stable histone states 
appear in a parameter regime whose size increases with increasing number of modification sites. We 
discuss the role of the spatial dependance, and we consider the effects of spatially heterogeneous 
enzymatic activity. Finally, we consider multistability in a model of several types of correlated 
post-translational modifications. 

PACS numbers: 



I. INTRODUCTION 

Nuclear chromosomes in eukaryotic organisms consist of the chromatin, a complex wrap that is primarily composed 
of DNA and histone proteins. The fundamental unit of the chromatin is the nucleosome, each of which contains 
two copies of the core histones H2A, H2B, H3 and H4, and approximately 150 base pairs of DNA. Each of the core 
histone proteins exhibits multiple amino acid residues that are subject to post-translational modifications (PTM) by 
chemical groups such as phospho-, acetyl-, methyl- or ubiquitin-groups that can be added and removed in a reversible 
manner. For example, H4 has a phosphorylation site, four acetylation sites and six methylation sites. Depending on 
the particular modification state of histones, certain regions of DNA in the chromatin are in an active or repressed 
state. Regulation of the PTMs of histones lies at the center of epigenetic control [IH3] . 

A commonly observed epigenetic phenomenon is the existence of alternative regulatory states. For example, in 
the fission yeast Schizosaccharomyces pombe the two mating type cassettes, mat2-P and mat3-M are usually in a 
silenced state in which the mating type genes are not expressed. When removing a portion of the silenced region and 
inserting a ura4+ reporter gene, the expression of ura4+ and the mating-type genes becomes bistable, with a state 
where ura4+ is repressed and a state where ura4+ is expressed [IHSj. The silenced state of ura4+ is associated with 
a high concentration of methylation marks on lysine of histone H3 (H3K9), while the active ura4+ state does not 
exhibit methylation of H3K9 7 . Each of the two epigenetic states is preserved under cell divisions, with transitions 
between them occuring only at a very low rate. 

Post-translational modifications arc regulated by various enzymes. In order to explain the appearance of multiple 
stable histone states, a non-local positive feedback mechanism has been put forward O [9] : A nucleosome that exhibits 
a particular modification recruits the enzymes that catalyze this modification. These enzymes then move to adjacent 
nucleosomes and cause the modification to be added there, a mechanism that has indeed been observed for some 
histone ace tyltransf erases, histone decacetylases and histone methyltransferases [101413) . Long-range feedback has been 
implemented in a stochastic simulation of a three-state model (unmodified state, acetylated state, methylated state) 
and it was shown to lead to robust bistability [TJ . Nearest-neighbour feedback has been considered in deterministic 
descriptions of two- and three-state models [111 HZ] • The authors of Ref. [TS] consider a two-state mean-field [T5] 
description that takes into account cooperativity in binding of enzymes, and they discuss the bifurcation diagram, 
including the effects of spatial dependence. In Ref. [17], the results of a stochastic simulation are compared to those of 
a mean-field description that does not explicitly consider spatial dependence. Perturbations due to cell divisions were 
considered, and instability of stable steady states due to such perturbations were found in the stochastic simulation, 
but not the mean-field approach. It is an open question how to obtain mean-field equations in the continuum 
starting from a stochastic description that predict the instabilities due to spatial dependance that are observed in the 
microscopic simulations. Among other things, this is one of the questions that we address in this work. 
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FIG. 1: Schematic illustration of an array of N nucleosomes, each of which contains S A = 5 PTMs A (blue) and S M = 4 
PTMs M (green) where PTMs of types A are regulated by a certain set of enzymes and PTMs M are regulated by another 
set of enzymes. Filled circles symbolize the presence of a PTM, empty circles indicate the absence of a PTM. In this example, 
occupations are = 5, n 1 ^ 1 = 0, = 1, n^ 1 = 2, etc. 



The considerable number of independently regulated modification sites in the chromatin has been hypothesized to 
give rise to a "histone code" [IS]: There are 2 T possible combinations of modified/unmodified configurations of T 
independently regulated PTMs, each of which potentially corresponds to a distinct "read-out" of information and 
ultimately a different epigenetic outcome. Recent efforts in identifying abundances of these histone modification 
states (also denoted as histone isoforms) have revealed that only few of the large number of possible isoforms are 
actually observed [TH1 HO]- It is also well known that regulation of different PTMs is correlated. For example, 
phosphorylation of H3 SerlO stimulates acetylation of H3 Lysl4 [31], and methylation of H3 Lys4 and Lys79 requires 
the ubiquitiniation of H2B Lysl23 [221 123] . In this work, we consider how such correlations in the regulation of PTMs 
reduce the information capacity of histone states. In particular, we study a model that is motivated by an interaction 
in the H3 N terminus where SerlO phosphorylation inhibits Lys9 methylation [24"] . 

We consider a master equation description of the stochastic dynamics of histone states (section [TTJ) . The system 
consists of a large number of nucleosomes, where each nucleosome exhibits several PTMs that are regulated by a par- 
ticular class of enzymes. We take into account the reversible addition and removal of PTMs due to enzymatic activity, 
as well as on-site ( "local" ) and nearest-neighbour ( "non-local" ) feedback mechanisms where modified nucleosomes re- 
cruit enzymes that either act locally or diffuse to adjacent nucleosomes. We use a quantum many-body formulation 
of the master equation a la Doi|25j and a tensor product state ansatz to obtain a system of nonlinear difference equa- 
tions (section III ) . We believe that the continuum limit of these equations is a suitable mean-field description that 
captures the role of spatial dependance in the master equation. The reader who is not interested in the derivation of 
the nonlinear difference equations/partial differential equations can go directly to Eqs. (10 1, Eqs. (13) and Eqs. (20 1. 
We numerically study the system of nonlinear partial differential equations (section |IV[). When considering one type 



of post-translational modification, and including at least two modification sites, bistable steady states are obtained 
without the necessity of explicit cooperativity at the level of the stochastic description (section [IV A ). The two stable 
steady states correspond to an unmodified state and a state with a high number of PTMs. We observe that increasing 
the number of modification sites increases the size of the parameter regime where bistable steady states exist. For 
a large number of modification sites, bistability is possible even if the coupling strength of the feedback mechanism 
is weak compared to the coupling strength of local processes. We observe that the spatial dependance due to the 
non-local feedback mechanism leads to instabilities of steady states under certain spatial perturbations of the histone 
state (section IV B I . These instabilities manifest themselves in traveling wave solutions of the system of nonlinear par- 
tial differential equations. We also consider spatially dependent rate parameters, which arise from adaptor proteins, 
such as DNA binding transcription factors, that recruit histone modifying enzymes to specific regions of chromatin 
(section IV C ) . We discuss how such spatially dependent enzyme activity gives rise to spatial heterogeneity in the 
epigenetic state. Finally, we introduce a model of two types PTMs that are regulated by different classes of enzymes 
and mutually inhibit each other (section IV D |. Such mechanisms are present in the chromatin, for example, in the 
case of H3 SerlO phosphorylation that inhibits H3 Lys9 methylation [23]. We find that inhibition in one direction is 
sufficient to reduce the full combinatorial set of four stable steady states to a set of three stable steady states where 
the presence of the two types of PTM is mutually exclusive. We conclude by discussing open problems and future 
directions. 



II. STOCHASTIC DYNAMICS OF HISTONE STATES 



We consider a one-dimensional array of N nucleosomes. Each nucleosome contains several modification sites of 
one or several independently regulated classes of PTMs, as schematically illustrated in Fig[T] A system comprised of 
N nucleosomes with S A modification sites of type A (e.g., acetylation) on each nucleosome is described by a state 
\n A , n A , rift) where the number of modified (e.g., acetylated) sites on nucleosome i is given by n A £ {0, 1, S A }. 
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We denote by P{n A ,n A , ...,nj^]t) the probability of finding the system in state \nf ,nf , ...,n^) at time t. In this 
and the following sections, we shall restrict ourselves to a single class of PTMs (i.e., regulated by a particular set of 
enzymes); however, in section [IV D| we shall discuss the case of two types of PTM. 

In the description of the stochastic dynamics of the histone state, we consider on-site ( "local" ) and nearest-neighbour 
("non-local" ) processes: 

1. The addition of a PTM A at nucleosome i with a rate X A , 

caused by enzymatic activity. 

2. The removal of a PTM A at nucleosome i with a rate /i A n A , 

A A 

A P 1 n i A -■ 

nf !-> nf - 1, 

as a result of enzymatic activity. 

3. The addition of a PTM A at nucleosome i with a rate f(nf_ 1 ,nf,nf +1 ), 

A ■<+!), Ail 

nf > nf + 1. 

The choice 

/(nf_!,nf,nf +1 ) = a A nf + a A (n A _ x + nf +1 - 2nf), (1) 

corresponds to a feedback mechanism that is both local and non-local. The first term (coupling parameter a A ) 
accounts for local feedback: the more PTMs are present at nucleosome i, the more enzymes that add PTMs of 
type A (e.g., acetylases) are present at «, and the more likely is the addition of further PTMs of type A. The 
second term (coupling parameter a A ) corresponds to non-local feedback: the enzymes at nearest-neighbouring 
nucleosomes i — 1 and i + 1 diffuse to nucleosome i and vice versa and, as in the case of local feedback, make 
the addition of additional PTMs more likely. 



4. The removal of a PTM A at nucleosome i with a rate nf (/(n^Lj^, nf , n A +1 ), i.e., 

A n?g(ntl, n ?, n ?+l) A -, 

n, > n; - 1. 



The choice 

g{nU, n f, n A +1 ) = p A {S A - n A ) + (3 A {2n A - n A _ x n A +1 ), (2) 

corresponds to a feedback mechanism that is both local and non-local. The first term (coupling parameter j3 A ) 
accounts for local feedback: The fewer PTMs are present at nucleosome i (i.e., the larger S — nf ), the more 
enzymes that cause the removal of PTM A (e.g., deacetylases) are present at i, making the removal of further 
PTMs more likely. The second term (coupling parameter (3 A ) corresponds to non-local feedback: the enzymes 
that cause the removal of PTMs A at nearest-neighbouring nucleosomes i — 1 and i + 1 diffuse to nucleosome i 
and vice versa and, as in the case of local feedback, make the removal of PTMs at site i more likely. 

The master equation for the above processes is given by 

dP(n A , n^; t) st^t^a . c i a a a i,, d / a a a , a a ,\ n , a a a a a .w 
-j. = 2^l A + /("«-!."! i«»+i)JL"(«i j ■■■,n i _ 1 ,n i - l,n i+ i, ...,n N ;t) - P(ni , ..,n,_i,n, ...,n N ;t)\ 

i=l 

N 

+ + g(nt-i, nf , nf +1 )][(nf + l)P(n A , ...,nf_ x ,nf + l,n i+1 , ...,nf,;t) - nfP(nf, nf-i, n t , nf + i, ...,n N ;t)]. (3) 
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III. DERIVATION OF NONLINEAR DIFFERENCE EQUATIONS 



We shall now introduce a notation of the master equation ((3| that is motivated by quantum physics [23 [26] . 
Standard quantum physics notation is used, i.e., \n A , ...,n^) = |nf ) ® |n^) (8 ... ® |^jv). We define 



A 4 



where the sum runs over all possible states. We introduce local raising and lowering operators [57] TZi and £j that 
are defined by 



K?\S,. 



0, 



; A iof) 



0, 



Indices A and i of operators signify that the operators are applied to state \n A ). When representing states |0 A ), 
|1 ),..., \S A ) by the S A + 1 unit vectors in S A + 1 dimensions, the lowering and raising operators can be represented 
by (S A + 1) x (S A + 1) dimensional matrices, 



/ 1 ... 
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The number operator is defined by M A = 1Z A £ A = Diag(0, 1, 2, S A ). In this notation, the master equation becomes 

0|*(t)) 



where % = H A <8> £ A ^ ... «. N t i, m 2 



9* 



where £•* denotes the S^-dimensional identity operator, and 



«!*(*)>, 

1 £^ (» . 



(4) 



T-Lfj (in simplified notation: "H = J2xLi^-t)^ 



Hf = \ A {n A -1 A ) + ^ A {C A -N A ) + (Tl A -l A )[ a A {N A l +N A 1 ~2N A ) + a A N A ' 
+(£? - K A )W A (Mti + M A - 2Mf) + fi A M% 



(5) 



where 1 A = Diag(l, 1, 1, 0), and M A = D\&g{S A ,S A - 1,...,1,0). In M, we substituted the functions |l| and 
([2]). We note that Q is an imaginary-time Schrodinger equation. The system corresponds to a quantum spin chain, 
though with a non-hermitian Hamitonian. 

The master equation Q is equivalent to a functional variation [28j . 



where 



6T 



r= / dt($\{d t -H)\V) 



(6) 



Since the system can be viewed as a quantum spin chain, albeit with a non-hermitian Hamiltonian "H, we make an 
ansatz for the wave-function in the Schrodinger picture as a tensor product state, 



N 



N 



i*(*)> =n !*<(*)>' 



(*i=n<*< 



(7) 



i=l 



and we write 1^,(4)) as a superposition of all possible states (we shall drop indices A from this point on), 



!*<(*)> = £cu*)l» 



n=0 



( CM \ 

Ci,i(*) 



V C i>s (t) J 



£< 

ra=0 



(8) 
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where X)n=o Q,n — 1j arL d Cj,n denotes the probability that nucleosome i has n modified sites. Since X)n=o = 1> 
this ansatz obeys the probabilistic constraint (^l^)^ n= o = 1 (i.e., expectation values ($|0|\&) of an observable O 
are properly normalized). 

Using this ansatz, the master equation in the formulation of (|6| becomes 



<9$ 



d *\ dCi >« /9 ^\H\*)) =0. (9) 



dCi. n / dt 



,fc=o 



Evaluating ^ yields a system of nonlinear difference equations for the probabilities C^ n that the nucleosome i has n 
modifications, 



dt 

dC i>n l<n<S 

dt 



-AG i)0 + fxC ia - G^clF? + &{m)) + + /3(m 4 )), 

-A(G ijn - Ci )n _i) - M(nCj,n - (n + l)C i)n+1 ) - (G i>n - C i) „_ 1 )(aF7 + a(n<)) 
-(nC<,„ - (n + l)G 4 ,„ +1 )(/3Gf + ^(m,)), 



= AC iiS -i-S/xC i)S + C i)S _ 1 (a^ v +a(n i ))-SC' i , s (^G7+ J 3(m i )), (10) 

where 

F? = {n,_i) - 2(m) + (n i+1 ) if 1< i < TV, F x v = -2<n x ) + (n 2 ) = (n w _ x ) - 2(n N ), 

Gf = K-x) - 2K) + <m i+1 ) if 1< i < AT, Gf 7 = -2<mi) + (ma) G^ = (mjv-i) - 2(mjv), 



and 



(n,) = ^nC,,„ (11) 

S 

(m.) = ^(5-n)G i , n = S-(n i ), (12) 



n=0 

(open boundary conditions). 

Equations ( 10 ) are a discretization of a system of nonlinear reaction-diffusion equations. Let £q be the lattice 
spacing (distance between nucleosomes) . In the mean-field/continuum limit a — > ch/£q, f3 — > (3/£q, and £o — » 0, we 
obtain the system of nonlinear partial differential equations for variables C n (x, t), n = 0, 1, S, 

^ = -AGo + ^-Go (aj^ (s^j + &X>Cr.)) + C * + -C«) - 

^ i<«<s _ x{Cn _ _ _ {n + 1)Cn+i) _ {(Jn _ Cn _j (a^f\ + a£>G s ) J 

\ s=0 ^ ' s=0 J 

-(nC n ~(n + l)C n+1 ) (^(iS-a)^) + jS£0S- *<?.)) ( 13 ) 

\ s=0 ^ X ' s=0 / 

^ = AC^x-^ + C^x f«E(^) +aj2(sC s ))-SC s Us - s)^) + ^(S - sC s )) . 

\ s=0 ^ ' ' s=0 / \ s=0 ^ ' s=0 / 

The diffusion terms are multiplied with the probabilities themselves. We note that the coefficient in front of the 
diffusion term is degenerate, and it is of interest to rigorously show the existence and stability of traveling wave 
solutions in reaction-diffusion equations of this type. 



IV. RESULTS 



In what follows, our analysis is based on numerical analysis of the system ( 10 ) over a finite parameter range. In the 
following, we set parameters a = 4a and j3 = 4/3, and we emphasize that varying the relative strength of local and 
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FIG. 2: Bifurcation diagram showing the steady state prob- 
abilities Co and Cs for S = 3 modification sites and param- 
eters A = ^t = l,/3 = 3asa function of parameter a. We 
denote the stable steady states where Co ~ 1 (few PTMs) 
and Cs ~ 1 (large number of PTMs) by X and Z, respec- 
tively, and the unstable steady state by Y. For a € [4.4, 7.5], 
steady states X, Y and Z appear, while for small a only X 
persists and for large a only Z persists. 



FIG. 3: Bifurcation diagram for S — 50 modification sites 
showing the steady state probabilities Ci ow = X!n=o Cn (i.e., 
low number of PTMs) and Chigh = 2n=46 ^™ (i- e -> high 
number of PTMs) of the stable steady states X and Z as a 
function of a. The remaining parameters are A = 5, fj, = 1, 
ft = 0.01. For a € [0.36, 0.53], bistability persists. Note that 
a <C A and ft <C fi. Inset: Width of the bistable regime in 
units of a as a function of the number of modification sites 
S (A = fj. = 1, /? = 3). It can be seen to increase linearly. 



non-local feedback does not qualitatively affect the results of our study. We note that as long as one is interested in 
the asymptotics (asymptotically long time) behavior of solutions of difference equations, what matters as input in the 
equations is the ratio (relative strength) of various coupling parameters (e.g., ft /a, A/a, etc.). One can always divide 
by a non-zero coupling parameters and rescale time to absorb this parameter in the left-hand-side of the difference 
equations. 

In section IV A, we will first discuss bistability in the model while neglecting spatial dependence. We will then 
incorporate spatial effects in part B, which we note fundamentally alters the picture. In section [IV C| we discuss the 
effects of spatial heterogeneity and in section [IV D| we discuss multiple correlated PTMs. 

A. Multiple stable steady states in the S-state model and the role of S 



In this section, we discuss the results of the nonlinear difference equations ( 10 ) when neglecting the spatial depen- 
dance, i.e., C\. n = Ci. n = ■■■ = Cm.u = C n . In this case, a system of coupled nonlinear ordinary differential equations 
(ODE) is obtained, 

dC 

= - AC + fid - 4aC 2 + iftCoCt , 
' • " 1 -= <S -A(C„ - C7„-i) - KnC n - (n + l)C n+l ) - AaC n (C n - C7„-i) - 40C n (nC n - (n + l)C„+i), 



dt 
dC s 
dt 



= XCs-i - SfiCs + AaCs-iCs - ASftC^. (14) 



Using this simplified ODE description, we evaluate steady states by setting dC n /dt = 0, and study their stability by 
analyzing the Jacobian matrix. Expressions for the steady state probabilities C n as a function of parameters A, fx, a 
and ft can be evaluated analytically. However, the resulting expressions are cumbersome and increasingly difficult to 
obtain for increasing S, and therefore calculations have been done numerically over a finite parameter range. 

For more than one modification site, i.e., S > 2, and appropriately chosen parameters (see below) we find that a 
parameter regime exists where three steady states coexist. The multistability is a consequence of the nonlinearities 
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in Eqs.(14) that are introduced by the feedback terms. Two of the steady states are stable attractors and one steady 
state is an unstable saddle point. We note that no explicit cooperativity is required in order to obtain bistability if S 
is chosen larger or equal than two. 

The bistability is illustrated in the bifurcation diagram of Fig. [2] where the steady state probabilities Co and C3 
are shown as a function of parameter a (the parameters used are S = 3, /x = A = l,/? = 3). If the feedback term 
for enzymes that catalyse the addition of PTMs is weak compared to the feedback term of enzymes that catalyse the 
removal of PTMs, only one steady state appears, as can be seen in Fig. [2] for a < 4.4. This steady state, which we 
denote by X, is characterized by Cq « 1, i.e., it corresponds to a state where very few PTMs are present. If the effects 
of the two terms that add PTMs approximately are roughly equal to the effects of the two terms that remove PTMs, 
three steady states exist (a £ [4.4, 7.5] in Fig. [2|. In addition to steady state X, a steady state with Cs ~ 1 appears. 
This steady state corresponds to a state with a high number of PTMs, and we shall denote it by Z. A third steady 
state (denoted by Y in Fig. [2]) is unstable. Finally, for large enough a, only steady state Z persists, as illustrated in 
Fig. D for a > 7.5. 

We note that in the previous paragraph we referred to the "strengths" of the four terms (1.-4. in section|n| as they 
can be read from the expectation values, e.g., ($| /(nj_i, rij, n^i)^}. In contrast, in the following paragraph, 
we shall refer to the magnitudes of the coupling parameters (i.e., a, /3, fj, 7 A) themselves. The values of the coupling 
parameters are controlled externally (e.g., the concentration, catalytic rate and diffusion rate of enzymes), while the 
expectation values also depend on system-dependent parameters (i.e., the number of modification sites S). 

Bistability is obtained only if both feedback terms are present, i.e., if both a and j5 are non-zero. If the number 
of modification sites, S, is small, bistable steady states appear only if the coupling parameters of the feedback terms 
are large compared to those of the local terms, i.e., only if the ratios A/a and /x//3 are small enough. However, 
with increasing number of modification sites S, the size of the parameter regime where multiple steady states appear 
increases, as shown in the inset of Fig. [3j and for large enough S, bistability can be established even if a <C A and 
/3 <C /i, as shown in Fig. [3] The existence of a large number of modification sites S that are regulated by a particular 
set of enzymes thus allows for a larger parameter regime of bistability. 



B. Spatial dependance 

In this section we will explicitly take into account spatial dependence, which is incorporated in the solutions to 



equations (10). We numerically integrate (10) and find that the stable steady states that were discussed in the 
previous section may become unstable for certain initial conditions. We illustrate this in Fig. [4] We set the initial 
probabilities Cj )W of the nucleosomes to those of steady state Z (the steady state where C^s is large), except for very 
few nucleosomes where we set the initial probabilities to values close to those corresponding to the second steady state 
X (35]. It can be seen that the system approaches steady state X, i.e., the spatially restricted perturbation of the 
histone state causes instability. This instability manifests itself by traveling wave solutions of the system of equations 
(10). It can be seen in Fig. [4] that for a perturbation away from the boundaries, two traveling wave fronts develop 
which travel at a constant velocity towards the boundaries of the system. If the perturbation is located at one of the 
boundaries of the system, only one wave front develops. 

There exists a set of parameters S, A, fi, a and /? where the velocity of the traveling wave(s) is zero. At that point, 
both steady states, X and Z, are stable with respect to spatial perturbations. For the parameters set of Fig. |4j this 
transition occurs at a* ~ 5.7 (bistability occurs for a G [4.4,7.5]). For a < a* and within range of bistability, the 
steady state X is the "stronger attractor" : If the initial state is Z and at least one nucleosome is perturbed such that 
its state is in the domain of fixed point X, the system approaches X, as is illustrated in Fig. [4] If the initial state is 
X, and at least one nucleosome is perturbed such that its state in the domain of steady state Z, the system bounces 
back into steady state X. In contrast, for a > a* , steady state Z is the "stronger attractor": If the initial state is X 
and at least one nucleosome is perturbed such that its state is in the domain of Z, the system approaches Z. If the 
initial state is Z, and at least one nucleosome is perturbed such that its state is in the domain of attraction of X, the 
system bounces back into steady state Z. 

In conclusion, for parameters a < a* , steady state X exhibits a very high degree of stability as any initial state 
of the system that gives rise to traveling wave solutions yields traveling waves that drive the system into state X. 
In contrast, for parameters a > a* , any traveling wave solution will drive the system into steady state Z. We note 



that when the asymptotic behaviour of equations ( 10 ) are considered, the number of nucleosomes in the system is not 
relevant. However, a larger number of nucleosomes does result in a longer duration for the traveling wave to spread 
over the entire system, which may be relevant if intermediate time scales are considered. 

Instabilitities due to traveling wave solutions could have significant impact on the stability and inheritance of 
chromatin steady states in daughter cells upon division. During cell division, it is thought that the parental nucle- 
osomes are randomly distributed among the two daughter cells, with the second half being newly synthesized [50] . 
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FIG. 4: Time evolution of probabilities d ,3 (t). The system 
(5 = 3 modification sites) is initially (time t = 0, red curve) 
in steady state Z (where C3 « 0.9), except for few nucleo- 
somes in the center that are strongly perturbed and whose 
probabilities are in the domain of steady state X. Parame- 
ters are \ = (i = 1, /3 = 3, a = 5.6. Two traveling wave 
fronts move towards the boundaries and drive the system 
into steady state X. The velocity of the waves is constant. 



FIG. 5: Time evolution of probabilities d,2,(t). Parameters 
are as in Fig. [4] except that A — 1 in part of the system, 
and A = 2 in the remainder, as indicated. The system is 
initially in steady state Z except for few nucleosomes in both 
regions whose states lie in the domain of X (red circles). In 
the left region (A = 1), two wave fronts move towards the 
boundaries, however, once the right front hits the A = 2 
region, it is stopped. In the A = 2 region, the perturbation 
does not cause the system to approach X. The reason is that 
for A = 1, steady state X is the "stronger attractor", while 
for A = 2, Z is the "stronger attractor" (terminology see 
text). 



The modification state of these new nucleosomes is crucial to the stability of the epigenetic state in the presence 
of non-local feedback terms. This can be seen as follows. The cell division can be modeled by replacing the states 
of half of the nucleosomes (randomly selected) at periodic intervals. Assume that the system is initially in steady 
state Z and parameters are set to the values of Fig. [4] where X is the "stronger attractor" . If the states of the newly 
synthesized nucleosomes are random (i.e., any state is possible), some of these nucleosomes might be in states that 
are in the domain of steady state X right after cell division. In this case, a traveling wave can form, and drive the 
system into steady state X (after one, several or many divisions, depending on the time-scales involved). We have 
verified this numerically. However, if the states of the newly synthesized nucleosomes are correlated with the state 
of the nucleosomes in the mother cell such that the states of the new nucleosomes are in the domain of attraction of 
the original state, such instabilities cannot arise. In the presence of non-local effects, a sufficient correlation between 
mother and daughter nucleosome states is hence necessary to preserve the chromatin state. This would relate to the 
notion of epigenetic memory and in fact there is a relation between daugher cell state and mother state [an example 
was discussed in the second paragraph of the introduction]. However, how this is conveyed at the molecular level 
remains a challenging open question. 

We conclude this section with a short discussion of the effects of considering explicit cooperative behaviour 
in the feedback terms. Explicit cooperative action of enzymes on-site, as well as of enzymes on nearest- 
neighbouring nucleosomes can be implemented using ansatz / coop (rti_i, rij, rij+i) = /(rij_i, Hi, Tii+i) + Srii—iriim+i 
and g coop (ni_i,ni, m+x) — g(rii-i, rii, n i+ i) +^f(S — m-i)(S — n,i)(S — n i+ i) Using the approach of sections [H III and 
|IV A| bistable steady states are observed, as was the case for the model without explicit cooperative action. However, 
bistability is possible even for the case S = 1. This in agreement with prior studies of two-state models with explicit 
cooperativity [121 [T7]. The difference equations that are obtained using this ansatz, or their continuum version, admit 
traveling wave solutions, as in the case of our model without explicit cooperative behaviour (10 1 where S > 2. 
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FIG. 6: Bifurcation diagram showing steady state probabilities C p , m of stable steady states for model (20 1 with two types of 
modifications, labeled by P and M, as a function of parameter /3 P_>M . Note that only P inhibits M, but not vice versa, i.e., 
pM^p = The remaining p aram eters are given by S p = S M = 2, A p = fi p = A M = fj, M = 1, j3 p = /3 M = 3, a p = a M = 4.5. 
For very small j3 p ^ M , four stable steady states (st.st.) exist, labeled by 00 (low P and low M), P0 (high P, low M), 0M (low 
P, high M), and PM (high P and high M). In an intermediate parameter regime, 

pP^,M e p.4 ! 3.3] i stable steady states 00, 
P0 and 0M persist, while for / g p ~ >M > 3.3 j only steady states 00 and P0 appear. We note that inhibition in only one direction 
(as p M ^ p = 0) is sufficient to obtain a parameter regime where steady states have either a high number of PTMs P or M, or 
neither, but not both. 



C. Spatially heterogeneous enzymatic activity 

In biological systems, nucleosome modifying enzymes are typically recruited to specific regions of the chromatin by 
adaptor proteins, such as DNA-binding transcription factors. As a result, the activity of these enzymes depends on 
the region of the chromatin. The increased or decreased activity of enzymes at certain nucleosomes can be taken into 
account by including a spatial dependance in parameters A and /i, i.e., A^ and /i^, where i is the nucleosome number. 
At each space point, the steady states are determined by the respective A^ and /Xj, i.e., the steady states locally 
correspond to the steady states with homogenous activity. Hence the parameter regimes where multiple stable steady 
states appear vary in size and position, and steady state probabilities Cj jW also depend on the nucleosome number 
i. For example, when choosing S = 3, fj, = 1, j3 = 3, and X = 1, bistability exists for a € [4.4,7.5] and a* « 5.7, 
while for parameters S — 3, /i = 1, /3 = 3 and A = 2, bistability persists for a £ [4.3,6.9], where a* « 5.5. As a 



consequence, for parameter a = 5.6, steady state X is the stronger attractor (in the sense explained in section IV B) 
for the former choice of parameters, while steady state Z is the stronger attractor for the latter choice of parameters. 
When perturbing a system that is initially in steady state Z in both A-regions, traveling wave solutions drive the 
system into steady state X at the nucleosomes where A = 1, but not in regions where A = 2, and the traveling waves 
in the region where A = 1 are stopped once they hit regions where A = 2, as shown in Fig. [5] Spatial dependence 
on the activity of histone modifying enzymes that is conferred by recruitment to regulatory regions of chromatin by 
transcription factors may thus stabilize the histone state from local and non-local perturbation. 



D. Multistability in model of several types of correlated PTMs 



Most proteins, such as histones, that are subject to PTM-dependent regulation are regulated via multiple modifica- 
tions. In this context, we discuss the results of including several types of modifications where each type is associated 
with different sets of enzymes, and thus different rate parameters A, fi, a, /3, a and f3. For example, one might 
consider different classes of acetylation (or phosphorylation, ubiquitination, etc.) sites, each of them associated with 
a different enzyme. Alternatively, one might consider PTMs of type P (e.g., phosphorylation) and PTMs of type M 
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(e.g., methylation) , with different rate parameters, A p and X M , a p and a M , etc. We denote by Ci tPim the probability 
of finding nucleosome i in the state with p PTMs of type P and m PTMs of type M. In this model, the number of 
stable steady states is four: the number of both M and P modifications is high (labeled by PM in the following), 
the number of P modifications is high and the number of M modifications is low (labeled by P0), the number of M 
modifications is high and the number of P modifications is low (labeled by 0M), and the number of both M and P 
modifications is low (labeled by 00). More generally, for T independent classes of modification sites, where a particular 
class of sites is associated with a particular set of coupling parameters, 2 T stable steady states are obtained. These 
steady states correspond to all possible combinations of states of high and low numbers of PTMs, i.e., all possible 
binary strings of length T. 

In practice, however, different types and sites of PTMs are often not independent from each other. There are 
examples where the presence of a certain PTM inhibits the addition of another PTM. An example is the H3 N 
tcminus where SerlO phosphorylation inhibits Lys9 methylation [24] . In the following, we derive difference equations 
using the formalism introduced in sections [TT] and |III| for a model of two types of PTMs, P and M, that mutually 
inhibit each other. We consider the processes 1.-4. (section [TTJ) separately for each of the two PTMs and add mutual 
inhibition (note that m = n M , p = n p ): 

n?nf nf(nf-l), (15) 

n P n M {n P _ 1)n M (16) 

In the case of ( 15 ), the presence of PTMs of type P leads to the removal of PTMs of type M, and in the case of ( |16[ ), 
the presence of PTMs of type M leads to the removal of PTMs of type P. 
The wave function is of form ^ with the local wave functions given by 

l*i(*)> =EE C itP , m (t)\p)\m) m = £ £ (plHe**», (17) 

p— m— p— m— 

where the normalization condition X) P =o Em=o ^i,p,m = 1 applies. The local operators Hi, Ci, Mi, Mi, li are defined 

as in section III and we denote the identity operator by £f~ (unity matrix of size S x ). Using this notation, the 
. . — „ „ . . . „ m - . 



"non-hermitian Hamiltonian" of the system is given by H — YIa—i Hi, where 

+Sf[X M + a M {M^ +K A ii ~ 2A^ M ) + a M M™\{Kf -T™) 
+[ M p + P P (MU + Mf +1 2M P ) + P P M P } (£f - N p )£^ 
+£ p [ M M + P M (M™ ! + Mf +1 - 2Mf ) + P M Mf](Cf - Nf 1 ) 

+/3 P^M (£ M _ tfMytfP + pM^P {c P _ tfPy/M. (18) 

The master equation in quantum variational formulation becomes 



5$ 



d* \dC itP , m I 0* =Q _ 



dt \ d<pi, p > 

Evaluating ( 19 ) yields a system of nonlinear difference equations for the probabilities Cj jPjm that the nucleosome at 
site i has p modifications of type P and m modifications of type M, 

= -(A p + a p F?r + a p (nf »(C i)P , m - C hP ^, n ) (X M + a M F?" + a M {nf ))(C itP , m C %p ^) 

-(pi p + (3 p Gj p + p p (m p ))( P C^ m - (p + l)C i>p+l , m ) - ( M M + (3 M Gj M + p M {mf ))(mC^ m - (m + l)C^ m+1 ) 
-P M ^ p (nf)(pC itP>m - (p+ l)C i>jH -i, m ) - /3 P ^ M (n M )(mC itPtm - (ro + l)C wn+1 ). (20) 

Here p = 0, 1, S p , m = 0, 1, S M , and 

F7 X ={nf- l )-Hnf) + {nf +l ) if 1< i < N, F? x = -2(raf ) + (n^) F** = K_ x ) - 2<n£), 

Gj x = (mf-x) - 2(mf ) + (m? +1 ) if 1 < i < N, Gj x = -2<mf ) + (m?) G V N * = (m*^) - 2<m$), 
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where X E {P,M}, and 



Sp=0 Em=0 pCi,p,m i 

s p S M 

Sp=0 Em=0 m Ci,p,m, , 



p=0 l~im=0\ J 



■p)Ci 
- m)Ci 



In Eqs.(20l, corrections for left-hand-side values of p = 0, 



= 0, f> = S, and m = S have to be taken into account, 
similarly as in the first and third equation of ([10 1. 

We consider the case of inhibition in only one direction by setting f3 M ^ p = and varying f3 p ^ M , as is the case 
in the example mention ed abo ve where SerlO phosphorylation inhibits Lys9 methylation. We evaluate steady states 
as explained in section (IVA|. For parameter choices of S p = S M — 2, A p = [i p = A M = fi M = 1, j3 p = /3 M = 3, 



a 



M 



4.5, a p = 4a p , a M = 4a M , f3 p = 4(3 P and /3 M = 4f3 M , we observe that for small fi p ~> M , all four stable 
steady states (as listed above) exist, as shown in Fig. [6} In an intermediate parameter regime only three stable steady 
states persist: 00, P0 and 0M, using the notation introduced above (Fig.pl). For large enough /3 P ^ M , only steady 
states 00 and P0 remain. This means that inhibitory interactions of two types of PTMs in only one direction are 
sufficient to obtain a parameter regime where steady states have either a high number of PTMs P or M, or neither, 
but not both. An analysis of traveling wave solutions of equations ( 20 ) similar to the one in section IV B applies in 
this case. 



V. CONCLUSIONS AND OUTLOOK 



The main results of this paper are as follows. We offer a robust method to obtain nonlinear partial differential 
equations describing the effective dynamics of histones. The method proceeds by mapping the system onto a quantum 
spin system whose dynamics is generated by a non-hermitian Hamiltonian. A feedback mechanism due to diffusion of 
enzymes along nucleosomes gives rise to multiple stable histone states. We study a number of novel aspects in histone 
systems that have not been reported before and are of biological relevance. We show that explicit cooperativity is 
not required to obtain multiple stable steady states as long as the number of PTMs is larger or equal to two, and we 
study the effects of varying the number of PTMs that are regulated by a particular set of enzymes. We also study 
the effect of spatially heterogeneous enzymatic on the histone state, and we apply our approach to a system of several 
correlated PTMs. 

Our approach can easily be generalized to higher spatial dimensions and more complicated network topologies. 
Processes other than the ones considered in this work could be included into the master equation and other biological 
systems might be studied. In the context of post-translational histone modifications, it might be of interest to consider 
more complex and more realistic systems. For example, the particular structure of the core histones might be taken 
into account i.e., the exact arrangement of the different modifications on the different core histones. Feedback processes 
among different types of post-translational modifications might be considered, as well as feedback loops that arise 
due to interactions between the histones and the DNA in the chromatin. It also remains an open question to study 
the existence and stability of traveling wave solutions in the nonlinear reaction-diffusion equations that arise in our 
model from a mathematically rigorous point of view. 

Acknowledgement. - We thank an anonymous referee for very helpful comments and suggestions that improved the 
presentation of the results in the paper. 
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